#######################################
# Summary Stats SHP
#######################################


rm(list = ls())
library(here)
library(data.table)
library(qs)
library(lfe)

library(tidyverse)
library(estimatr)
library(texreg)
library(ggthemes)
library(arsenal)
library(haven)

source(here("code", "functions.R"))
shp <- qread(here("data", "shp.qs")) 



# Controls ----------------------------------------------------------------
shp <- shp %>% 
  mutate(
    income = replace(hhinc_net, hhinc_net <0 , NA),
    income = income/12,
    income = case_when(income < 3000 ~ "$<$ CHF3000", 
                       income <= 5000 ~ "CHF3000-5000", 
                       income <= 7000 ~ "CHF5000-7000", 
                       income <= 9000 ~ "CHF7000-9000",
                       income > 9000 ~ "$>$ CHF9000"),
    income = factor(income, levels = c("$<$ CHF3000",
                                       "CHF3000-5000",
                                       "CHF5000-7000",
                                       "CHF7000-9000",
                                       "$>$ CHF9000")),
    age = as.numeric(age),
    female = recode(as.numeric(gender), `1` = 0, `2` = 1, .default = NA_real_),
    education = recode(as.numeric(educat), `0` = "Less than primary", `1` = "Primary",
                       `2` = "Secondary", `3` = "Secondary", `4` = "Secondary",
                       `5` = "Secondary", `6` = "Secondary", `7` = "Tertiary",
                       `8` = "Tertiary", `9` = "Tertiary", `10` = "Tertiary", .default = NA_character_),
    education = factor(education, levels = c("Less than primary", "Primary", "Secondary", "Tertiary")),
    employment = recode(as.numeric(wstat), `1` = "Employed", `2` = "Unemployed", `3` = "Not in labor force", .default = NA_character_),
    employment = factor(employment, levels = c("Employed", "Unemployed", "Not in labor force")),
    married = recode(as.numeric(civsta), `1` = 0, `2` = 1, `3` = 0, `4` = 0, `5` = 0, `6` = 0, `7` = 0, .default = NA_real_),
    opin_foreign = recode(as.numeric(opin_foreign), `1` = "Favor equality with foreigners", `2` = "Neither", `3` = "Favor opportunities for Swiss", .default = NA_character_),
    opin_foreign = factor(opin_foreign, levels = c("Favor equality with foreigners", "Neither", "Favor opportunities for Swiss")),
    polit_interest = case_when(as.numeric(polit_interest) >= 0 ~ as.numeric(polit_interest), TRUE ~ NA_real_),
    languagereg = recode(languagereg2019, `1` = "German", `2` = "French", `3` = "Italian", `4` = "Romansh")
)

#Filter to main sample
df <- shp %>% filter(BR == 1 & travelMUNmin <= 30 & complete.cases(weight_i_pop))

#Define treatmetn and control
df <- df %>% 
  mutate(treated = as.numeric(travelMUNmin <= 15),
         treated = recode(treated, `1` = "Treated (0-15 min)", `0` = "Control (15-30 min)"),
         treated = factor(treated)) %>% 
  mutate(weights = weight_i_pop * n() / sum(df$weight_i_pop))

# Descriptive Stats 1999 -------------------------------------------------------
shp_summary1999 <- tableby(treated ~ income + age + female + education + employment + married + opin_foreign + polit_interest + languagereg, 
                       weights = weight_i_pop, data = df %>% filter(complete.cases(weight_i_pop), year == 1999),
                      digits = 1, numeric.stats =  c("Nmiss2", "meansd", "range")) %>% 
  summary(text = "latex",
          labelTranslations = list(income = "Household Income",
                                   age = "Age (Years)",
                                   female = "Female",
                                   education = "Highest Education Level",
                                   employment = "Employment",
                                   married = "Currently Married",
                                   opin_foreign = "Opinion on Foreigners",
                                   polit_interest = "Political Interest",
                                   languagereg = "Language Region"),
          pfootnote = T, width = 5, fontsize = "scriptsize")

shp_summary1999

capture.output(shp_summary1999, file= here("tables", "summary_stats", "shp_summary1999.tex"))



# Descriptive Stats 2017 -------------------------------------------------------
shp_summary2017 <- tableby(treated ~ income + age + female + education + employment + married + opin_foreign + polit_interest + languagereg, 
                           weights = weight_i_pop, data = df %>% filter(complete.cases(weight_i_pop), year == 2017),
                           digits = 1, numeric.stats =  c("Nmiss2", "meansd", "range")) %>% 
  summary(text = "latex",
          labelTranslations = list(income = "Household Income",
                                   age = "Age (Years)",
                                   female = "Female",
                                   education = "Highest Education Level",
                                   employment = "Employment",
                                   married = "Currently Married",
                                   opin_foreign = "Opinion on Foreigners",
                                   polit_interest = "Political Interest",
                                   languagereg = "Language Region"),
          pfootnote = T, width = 5, fontsize = "scriptsize")

capture.output(shp_summary2017, file= here("tables", "summary_stats", "shp_summary2017.tex"))




# Descriptive Stats 1999 -------------------------------------------------------
shp_summary_TI1999 <- tableby(treated ~ income + age + female + education + employment + married + opin_foreign + polit_interest, 
                           weights = weight_i_pop, data = df %>% filter(complete.cases(weight_i_pop), year == 1999, canton_name == "TI"),
                           digits = 1, numeric.stats =  c("Nmiss2", "meansd", "range")) %>% 
  summary(text = "latex",
          labelTranslations = list(income = "Household Income",
                                   age = "Age (Years)",
                                   female = "Female",
                                   education = "Highest Education Level",
                                   employment = "Employment",
                                   married = "Currently Married",
                                   opin_foreign = "Opinion on Foreigners",
                                   polit_interest = "Political Interest",
                                   languagereg = "Language Region"),
          pfootnote = T, width = 5, fontsize = "scriptsize")

shp_summary_TI1999

capture.output(shp_summary_TI1999, file= here("tables", "summary_stats", "shp_summary_TI1999.tex"))



# Descriptive Stats 2017 -------------------------------------------------------
shp$polit_interest[shp$year == 2017] %>% summary
shp_summary_TI2017 <- tableby(treated ~ income + age + female + education + employment + married + opin_foreign + polit_interest, 
                           weights = weight_i_pop, data = df %>% filter(complete.cases(weight_i_pop), year == 2017, canton_name == "TI"),
                           digits = 1, numeric.stats =  c("Nmiss2", "meansd", "range")) %>% 
  summary(text = "latex",
          labelTranslations = list(income = "Household Income",
                                   age = "Age (Years)",
                                   female = "Female",
                                   education = "Highest Education Level",
                                   employment = "Employment",
                                   married = "Currently Married",
                                   opin_foreign = "Opinion on Foreigners",
                                   polit_interest = "Political Interest",
                                   languagereg = "Language Region"),
          pfootnote = T, width = 5, fontsize = "scriptsize")

capture.output(shp_summary_TI2017, file= here("tables", "summary_stats", "shp_summary_TI2017.tex"))

